## load packages
library(lmerTest)
library(lmtest)
library(lme4)
library(car)
library(effects)
library(gplots)
library(plotrix)
library(psych)
library(ggplot2)
library(grid)
library(gridExtra)
library(scales)
library(MASS)
library(lsmeans)
library(pbkrtest)
library(cowplot)
library(devtools)
library(dplyr)
#install_github("vqv/ggbiplot")
#library(ggbiplot)

Background

Data is from Kāne’ohe Bay, O’ahu, Hawai’i. Data extends across two archipelagic bleaching events and the subsequent recovery periods where tempertaure stress subsided. Data periods represent:
- first bleaching (October 2014) - first recovery (February 2015) - second bleaching (October 2015)
- second recovery (February 2016)

Physiological and immunological parameters were collected from from Montipora capitata collected from two locations (Sites):
- Lilipuna – a shallow fringing reef west of the Hawaii Insititute of Marine Biology
- Reef 14 – a shallow inshore patch reef in central Kāne’ohe Bay

Lilipuna has been categorized as experiencing near shore impacts from freshwater inundation, subteranean groundwater discharge, and seawater biogeochemistry influeces by long residence times, low flow, and riverine/terrigenous nutrient inputs. Additionally, pCO2 flux at this site is “moderate”, with on average ~200ppm pCO2 +/- flux

Reef 14 has seawater with a shorter residence time, little terrestrial impact, but has a greater pCO2 flux due to reef metabolism. This flux can be substantial, being > 600 ppm pCO2 in a 24h period and reaching maximum values of >900ppm pCO2.

Import and Normalize

Raw data from physiological assessments is imported and normalized according to established protocols. Immunolgical data has been normalized previously and is represented as the final dataframe.

###########################################################################
## data period: October 2014, February 2015, October 2015, February 2016 ##
###########################################################################

# working directory
main<-setwd(getwd())
rm(list=ls())

### data file
# all lab data from all time points (physio and immuno, PLUS color analysis)
physimmun<-read.csv("data/Gates_Mydlarz_20142016_physimmun.csv", header=T, na.string=NA) #physimmuno data
qPCR.data<-read.csv("data/Gates_Mydlarz_20142016_qPCR.csv", header=T, na.string=NA) #qPCR data

data<-read.csv("data/Gates_Mydlarz_20142016_ALL_DATA.csv", header=T, na.string=NA) #qPCR data


################################################
# PHYSIOLOGY: calculate, normalize dependent variables
################################################
data$cells..ml<-as.numeric(data$cells..ml)
data$ID<-as.factor(data$ID)


# helpful shorthand
SA<-data$surface.area.cm2 # surface area in cm2
blastate<-data$total.blastate.ml # tissue slurry blastate in ml

# Symbiodinium.cells..cm2 == cell.ml * blastate / SA
data$zoox<- (data$cells..ml*blastate)/SA

# ug.chlorophyll.a..cm2 == ug.chl.a.ml * blastate / SA
data$chla<- (data$ug.chla..ml*blastate)/SA

# pg.chlorophyll.a..cell == ug.chl.a.ml * 10^6 / cells.ml
data$chlacell<- (data$ug.chla..ml*10^6)/data$cells..ml

# AFDW.mg..cm2 == convert AFDW g to mg, mutiply by blastate volume, divide by cm2
data$AFDW<- (data$g.AFDW..ml*1000*blastate)/SA

# mg.protein..cm2 == mg.protein.ml * blastate / SA
data$prot<- (data$mg.prot..ml*blastate)/SA


#### rearrange and clean up

# drop unwanted columns by name
data<-data[!colnames(data) %in% c("total.blastate.ml", "surface.area.cm2", "mg.prot..ml", "cells..ml", "ug.chla..ml", "g.AFDW..ml")]

data<-data[!(data$ID=="N37"), ] # this coral has an NA for symbiont type, remove from df

# reorder columns to show: hierarchy of structure, physio., immuno., colorimetry
# this is now the working dataframe
df<-data[, c(1:6, 20:24, 7:19)]

# remove negative values for MEL and POX and
df$MEL[df$MEL <= 0]<- NA
df$MEL[df$MEL > 0.05]<- NA

df$POX[df$POX <= 0]<- NA
df$CAT[df$CAT <= 0]<- NA
df$PPO[df$PPO <= 0]<- NA


df$chlacell[df$chlacell >= 10]<- NA # removing 2 outliers
df$CAT[df$CAT >= 900]<- NA # 2 outliers removed


# don't need "Ramp" lab data, this was lab thermally stressed. Remove this.
df<-df[!(df$Status=="Lab-Ramp"),]

#write.csv(df, "df.normalized.csv")

Color scores

Overall, color is a poor measure of bleaching in these fragments, and there is not a strong relationship between color (R/G/B) and physiological bleaching quantification. Issues may be in the approach to get color scores (in shooting images, or downstream in photoshop), as well as in the fact the same color score can represent significant changes in cell density or chla content thereby making the relationship noisy (although significant).
- red and green are best explantory variable for chla (~25 % R2)
- blue and green slightly better to explain zoox density (~33 % R2)

Color score by Chlorophyll a cm-2

par(mfrow=c(1,3), mar=c(5,4,1,2), pty="sq")

##########################################
# testing color scores and chla relationship
mod<-lm(coral.red.adj~chla, data=df)
plot(coral.red.adj~chla, data=df)
abline(mod, col="red")

mod<-lm(coral.blue.adj~chla, data=df)
plot(coral.blue.adj~chla, data=df)
abline(mod, col="blue")

mod<-lm(coral.green.adj~chla, data=df)
plot(coral.green.adj~chla, data=df)
abline(mod, col="green")


### Color score by *Symbiodinium* cells cm-2
##########################################
# testing color score and zoox relationship
mod<-lm(coral.red.adj~zoox, data=df)
plot(coral.red.adj~zoox, data=df)
abline(mod, col="red")

mod<-lm(coral.blue.adj~zoox, data=df)
plot(coral.blue.adj~zoox, data=df)
abline(mod, col="blue")

mod<-lm(coral.green.adj~zoox, data=df)
plot(coral.green.adj~zoox, data=df)
abline(mod, col="green")


####### PCA and color score
####### Run the PCA first and save PC1 data 

pca.df<-(df[-c(1:113),c(2:6, 17:19)]) # all data with only "Event" and "Site" as categories
colnames(pca.df)<-c("Period", "Status", "Site", "Status_Site", "ID", "red", "green", "blue")
pca.df<-pca.df[-9, ] #remove NA row

# apply PCA - scale. = TRUE for center as advised
PC <- prcomp(pca.df[, 6:8], center = TRUE, scale.= TRUE)  

#correlatin matrix helps to standardize the data and account for different scales
# print(PC) # to print PC loadings

# summary(PC) #signficance of PCs: proportion of variance captured, cumulative variance
# plot(PC, type="lines") will plot PCs

newdat<-PC$x[,1:3]; # newdata frame with only PCs 
# PC1 explain ~>82% of variance

# save the file of PCs
write.csv(newdat, "PC1-allcolors_alltimes.csv")

# the SDEV^2 is the eignevalue
ev<-PC$sdev^2 # PC1-3 have eignevalues >1
# ev/sum(ev) # to give proportional eignevalues

#########
# plot it as PCA analysis by bleaching period

## PC1 and PC2
g <- ggbiplot(PC, choices = 1:2, obs.scale = 1, var.scale = 1, 
              groups= pca.df[,1], ellipse = TRUE, 
              circle = FALSE) +
              scale_color_discrete(name = '') +
  theme_bw() +
  theme(legend.text=element_text(size=15)) +
  theme(panel.background = element_rect(colour = "black", size=1))+
  theme(legend.key = element_blank())+
  theme(legend.direction = 'horizontal', legend.position = 'top') +theme(aspect.ratio=0.7)
print(g)

##################

### graphing the PC1 color score relationship with chlorophyll or *Symbiodinium* density**
### significant relationships, but only 27% (chla) and 35% (zoox) of variance explained

##########################################
# testing color scores and chla relationship
mod<-lm(PC1~chla, data=df)
plot(PC1~chla, data=df)
abline(mod, col="mediumorchid3")

##########################################
# testing PC1 score and zoox relationship
mod<-lm(PC1~zoox, data=df)
plot(PC1~zoox, data=df)
abline(mod, col="mediumorchid")



### Categorical binning for bleaching
# First, looked at all the data as histograms of chlorophyll a and zoox density across all 4 time periods. Second, looked at the histograms forjust bleaching, and then just recovery periods. While it is difficult to say when a coral is truly bleached, a conservative measure of < 3 ug chla/cm2 or < 1.5 million cells may be used as a relative measure of bleaching.


#############################
# all data: bleaching and recovery periods
par(mfrow=c(1,2), mar=c(5,4,1,2), pty="sq")
hist(df$zoox, main="all time points"); hist(df$chla, main="all time points")


#############################
#### just "bleaching" periods
bleaching.df<-df[df$Status=="Bleaching", ]

par(mfrow=c(1,2), mar=c(5,4,1,2), pty="sq")
hist(bleaching.df$zoox, main="bleaching periods")
hist(bleaching.df$chla, main="bleaching periods") 
# less than 3 ug chla/cm2 is good indication of "bleached"
# less than 1.5 million cells/cm2 is good indication of "bleached

# what about symbiont community
hist(bleaching.df[bleaching.df$dom=="D",]$chla)
hist(bleaching.df[bleaching.df$dom=="C",]$chla)


#############################
#### just "recovery" periods
recovery.df<-df[df$Status=="Recovery", ]

par(mfrow=c(1,2), mar=c(5,4,1,2), pty="sq")
hist(recovery.df$zoox, main="recovery periods")
hist(recovery.df$chla,  main="recovery periods")

### Bleaching categories
# Decided to bin the data based on bleaching (< 40%  ug chla cm-2) and pigmented (> 60% ug chla cm-2). This should help in finding those corals more affected by the thermal stress, and those not affected as severely.

#### applying binning

#### no break here for Lab 2014--these are all effectively field controls for physiology
Field.2014<-df[(df$Period=="2014 Field.Feb"),]
Field.2014<- Field.2014 %>% mutate(
            bleach.category = "Field-Pre")

#### no break here for Lab 2014--these are all effectively field controls for physiology
Lab.2014<-df[(df$Period=="2014 Lab"),]
quantile(Lab.2014$chla, 0.4, na.rm=TRUE) # = 3.86
Lab.2014<- Lab.2014 %>% mutate(
            bleach.category = "Lab-pigmented")

#### 40% quantile for October 2014
Oct.2014<-df[(df$Period=="2014 Oct"),]
quantile(Oct.2014$chla, 0.4, na.rm=TRUE) # = 2.97
Oct.2014<- Oct.2014 %>% mutate(
            bleach.category = ifelse(chla < "2.97", "pale", "pigmented"))

#### 40% quantile for February 2015
Feb.2015<-df[(df$Period=="2015 Feb"),]
quantile(Feb.2015$chla, 0.4, na.rm=TRUE) # = 3.72
Feb.2015<- Feb.2015 %>% mutate(
            bleach.category = ifelse(chla < "3.72", "pale", "pigmented"))

#### 40% quantile for October 2015
Oct.2015<-df[(df$Period=="2015 Oct"),]
quantile(Oct.2015$chla, 0.4, na.rm=TRUE)  # = 2.49
Oct.2015<- Oct.2015 %>% mutate(
            bleach.category = ifelse(chla < "2.49", "pale", "pigmented"))

#### 40% quantile for January 2016
Feb.2016<-df[(df$Period=="2016 Feb"),]
quantile(Feb.2016$chla, 0.4, na.rm=TRUE)  # = 3.84
Feb.2016<- Feb.2016 %>% mutate(
            bleach.category = ifelse(chla < "3.84", "pale", "pigmented"))


##########
## Now combine all the dataframes above back into a single df
df<-rbind(Field.2014, Lab.2014, Oct.2014, Feb.2015, Oct.2015, Feb.2016)

Summary dataframes

These dataframes show mean, SE and n for each group. This can be subset to make figures or exported as a summary file. The header of the dataframe is shown below…

# make summary dataframes

# dataframes of means can later be divided into categories by "C or D" dominated

###------#### means

# physiology
zoox.m<-aggregate(zoox~Site+Status+Period+dom,data=df, mean)
chla.m<-aggregate(chla~Site+Status+Period+dom,data=df, mean)
chlacell.m<-aggregate(chlacell~Site+Status+Period+dom, data=df, mean)
prot.m<-aggregate(prot~Site+Status+Period+dom, data=df, mean)
AFDW.m<-aggregate(AFDW~Site+Status+Period+dom, data=df, mean)

# imunology
CAT.m<-aggregate(CAT~Site+Status+Period+dom,data=df, mean)
POX.m<-aggregate(POX~Site+Status+Period+dom,data=df, mean)
SOD.m<-aggregate(SOD~Site+Status+Period+dom,data=df, mean)
PPO.m<-aggregate(PPO~Site+Status+Period+dom,data=df, mean)
MEL.m<-aggregate(MEL~Site+Status+Period+dom,data=df, mean)

###------#### SE
# physiology
zoox.SE<-aggregate(zoox~Site+Status+Period+dom,data=df, std.error)
chla.SE<-aggregate(chla~Site+Status+Period+dom,data=df, std.error)
chlacell.SE<-aggregate(chlacell~Site+Status+Period+dom, data=df, std.error)
prot.SE<-aggregate(prot~Site+Status+Period+dom, data=df, std.error)
AFDW.SE<-aggregate(AFDW~Site+Status+Period+dom, data=df, std.error)

# imunology
CAT.SE<-aggregate(CAT~Site+Status+Period+dom,data=df, std.error)
POX.SE<-aggregate(POX~Site+Status+Period+dom,data=df, std.error)
SOD.SE<-aggregate(SOD~Site+Status+Period+dom,data=df, std.error)
PPO.SE<-aggregate(PPO~Site+Status+Period+dom,data=df, std.error)
MEL.SE<-aggregate(MEL~Site+Status+Period+dom,data=df, std.error)

###------#### n-value
# physiology
zoox.n<-aggregate(zoox~Site+Status+Period+dom,data=df, length)
chla.n<-aggregate(chla~Site+Status+Period+dom,data=df, length)
chlacell.n<-aggregate(chlacell~Site+Status+Period+dom, data=df, length)
prot.n<-aggregate(prot~Site+Status+Period+dom, data=df, length)
AFDW.n<-aggregate(AFDW~Site+Status+Period+dom, data=df, length)

# imunology
CAT.n<-aggregate(CAT~Site+Status+Period+dom,data=df, length)
POX.n<-aggregate(POX~Site+Status+Period+dom,data=df, length)
SOD.n<-aggregate(SOD~Site+Status+Period+dom,data=df, length)
PPO.n<-aggregate(PPO~Site+Status+Period+dom,data=df, length)
MEL.n<-aggregate(MEL~Site+Status+Period+dom,data=df, length)


#### physio dataframe
fig.df.phys<-cbind(zoox.m, zoox.SE[c(5,0)],chla.m[c(5,0)], chla.SE[c(5,0)], chlacell.m[c(5,0)], chlacell.SE[c(5,0)], prot.m[c(5,0)], prot.SE[c(5,0)], AFDW.m[c(5,0)], AFDW.SE[c(5,0)], chla.n[c(5,0)])

colnames(fig.df.phys)<-c("Site","Status", "Period", "dom", "zoox", "zoox.SE", "chla", "chla.SE", "chlacell", "chlacell.SE", "prot", "prot.SE", "AFDW", "AFDW.SE", "N-chla")

# make an interaction column
fig.df.phys$int<-interaction(fig.df.phys$Site, fig.df.phys$dom)
fig.df.phys$int<-factor(fig.df.phys$int)
####


#### immuno dataframe
fig.df.immuno<-cbind(CAT.m, CAT.SE[c(5,0)], POX.m[c(5,0)], POX.SE[c(5,0)], SOD.m[c(5,0)], SOD.SE[c(5,0)], PPO.m[c(5,0)], PPO.SE[c(5,0)], MEL.m[c(5,0)], MEL.SE[c(5,0)])

colnames(fig.df.immuno)<-c("Site","Status", "Period", "dom", "CAT", "CAT.SE", "POX", "POX.SE", "SOD", "SOD.SE", "PPO", "PPO.SE", "MEL", "MEL.SE")

# make an interaction column
fig.df.immuno$int<-interaction(fig.df.immuno$Site, fig.df.immuno$dom)
fig.df.immuno$int<-factor(fig.df.immuno$int)

# write.csv(xxx, "Figure.df.xxx.csv")

Figures

# side by side figures separated by Sites 
### Physiological figures

color.scheme<-c("gray80", "gray60", 'lightskyblue', 'dodgerblue', 
                "gray50", "gray30", "thistle", "coral")
break.points<-c("Lilipuna.Pre.C", "Lilipuna.Pre.D", "Lilipuna.C", "Lilipuna.D", 
                "Reef 14.Pre.C", "Reef 14.Pre.D", "Reef 14.C", "Reef 14.D")
legend.names<-c("Lil.Pre C", "Lil.Pre D", "Lil C", "Lil D", 
                "R14.Pre C", "R14.Pre D", "R14 C", "R14 D")
axis.labels<-c("Feb 2014 \nPre-Bleaching", "Oct 2014 \nBleaching","Feb 2015 \nRecovery","Oct 2015 \nBleaching", "Feb 2016 \nRecovery")


Fig.formatting<-(theme_classic()) +
  theme(text=element_text(size=10),
  axis.line=element_blank(),
  legend.text.align = 0,
  legend.text=element_text(size=10),
    #legend.title = element_blank(),
      panel.border = element_rect(fill=NA, colour = "black", size=1),
      aspect.ratio=1, 
    axis.ticks.length=unit(0.25, "cm"),
    axis.text.y=element_text(
      margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=10), 
    axis.text.x=element_text(
      margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=8)) +
  theme(legend.key.size = unit(0.4, "cm")) +
  theme(aspect.ratio=1) +
  theme(panel.spacing=unit(c(0, 0, 0, 0), "cm"))


pd <- position_dodge(0.15) #offset for error bars


########################
#### graph as category of C/D
######################## 

fig.df.phys$int<-factor(fig.df.phys$int, levels=c("Lilipuna.Pre.C", "Lilipuna.Pre.D", "Lilipuna.C", "Lilipuna.D", "Reef 14.Pre.C", "Reef 14.Pre.D", "Reef 14.C", "Reef 14.D"))

###################
# zoox 
zoox.fig<-ggplot(data=fig.df.phys, aes(x=Period, y=zoox/10^6, group=int, color=int)) + geom_errorbar(aes(ymin=zoox/10^6-zoox.SE/10^6, ymax=zoox/10^6+zoox.SE/10^6), size=.8, width=0, position=pd) +
  geom_line(position=pd, size=.8) +
  geom_point(position=pd, size=3, pch=19) +
  coord_cartesian(ylim=c(0, 3)) +
  ylab(expression(paste(~italic("Symbiodinium") ~(10^6~cells~cm^-2), sep=""))) +
  scale_color_manual(values=c(color.scheme),
                     name= "Site:Clade",
                     breaks=break.points,
                     labels=legend.names) +
  scale_x_discrete(name ="Year and Event", 
      labels=axis.labels)+
  Fig.formatting

zoox.fig

Chlorophyll a cm-2

########## chla cm2 ############
chla.fig<-ggplot(data=fig.df.phys, aes(x=Period, y=chla, group=int, color=int)) + geom_errorbar(aes(ymin=chla-chla.SE, ymax=chla+chla.SE), size=.8, width=0, position=pd) +
  geom_line(position=pd, size=.8) +
  geom_point(position=pd, size=3) +
  coord_cartesian(ylim=c(0, 8)) +
  ylab(expression(paste("Chlorophyll", ~italic("a"), ~(mu*g~cm^-2), sep=""))) +
  scale_color_manual(values=c(color.scheme),
                     name= "Site:Clade",
                     breaks=break.points,
                     labels=legend.names) +
  scale_x_discrete(name ="Year and Event", 
      labels=axis.labels)+
  Fig.formatting
chla.fig

Chlorophyll a cell-2

########## chla cell ############
chlacell.fig<-ggplot(data=fig.df.phys, aes(x=Period, y=chlacell, group=int, color=int)) + geom_errorbar(aes(ymin=chlacell-chlacell.SE, ymax=chlacell+chlacell.SE), size=.8, width=0, position=pd) +
  geom_line(position=pd, size=.8) +
  geom_point(position=pd, size=3) +
  coord_cartesian(ylim=c(0, 8)) +
  ylab(expression(paste("Chlorophyll", ~italic("a"), ~(pg~cell^-1), sep=""))) +
  scale_color_manual(values=c(color.scheme),
                     name= "Site:Clade",
                     breaks=break.points,
                     labels=legend.names) +
  scale_x_discrete(name ="Year and Event", 
      labels=axis.labels)+
  Fig.formatting
chlacell.fig

Protein biomass cm-2

########## protein ############
prot.fig<-ggplot(data=fig.df.phys, aes(x=Period, y=prot, group=int, color=int)) + geom_errorbar(aes(ymin=prot-prot.SE, ymax=prot+prot.SE), size=.8, width=0, position=pd) +
  geom_line(position=pd, size=.8) +
  geom_point(position=pd, size=3) +
  coord_cartesian(ylim=c(0, 1)) +
  ylab(expression(paste("Protein", ~(mg~cm^-2), sep=""))) +
  scale_color_manual(values=c(color.scheme),
                     name= "Site:Clade",
                     breaks=break.points,
                     labels=legend.names) +
  scale_x_discrete(name ="Year and Event", 
      labels=axis.labels)+
  Fig.formatting
prot.fig

Total biomass (AFDW) cm-2

########## AFDW ############
AFDW.fig<-ggplot(data=fig.df.phys, aes(x=Period, y=AFDW, group=int, color=int)) + geom_errorbar(aes(ymin=AFDW-AFDW.SE, ymax=AFDW+AFDW.SE), size=.8, width=0, position=pd) +
  geom_line(position=pd, size=.8) +
  geom_point(position=pd, size=3) +
  coord_cartesian(ylim=c(0, 60)) +
  ylab(expression(paste("Total biomass", ~(mg~cm^-2), sep=""))) +
  scale_color_manual(values=c(color.scheme),
                     name= "Site:Clade",
                     breaks=break.points,
                     labels=legend.names) +
  scale_x_discrete(name ="Year and Event", 
      labels=axis.labels)+
  Fig.formatting
AFDW.fig

Imunological figures

Prophenoloxidase (PPO)

#names(fig.df.immuno)

# site means
CAT.m2<-aggregate(CAT~Site+Status+Period,data=df, mean)
POX.m2<-aggregate(POX~Site+Status+Period,data=df, mean)
SOD.m2<-aggregate(SOD~Site+Status+Period,data=df, mean)
PPO.m2<-aggregate(PPO~Site+Status+Period,data=df, mean)
MEL.m2<-aggregate(MEL~Site+Status+Period,data=df, mean)

###------#### SE
# imunology
CAT.SE2<-aggregate(CAT~Site+Status+Period,data=df, std.error)
POX.SE2<-aggregate(POX~Site+Status+Period,data=df, std.error)
SOD.SE2<-aggregate(SOD~Site+Status+Period,data=df, std.error)
PPO.SE2<-aggregate(PPO~Site+Status+Period,data=df, std.error)
MEL.SE2<-aggregate(MEL~Site+Status+Period,data=df, std.error)

immuno.site.df<-cbind(CAT.m2, CAT.SE2[c(4,0)], POX.m2[c(4,0)], POX.SE2[c(4,0)], SOD.m2[c(4,0)], SOD.SE2[c(4,0)], PPO.m2[c(4,0)], PPO.SE2[c(4,0)], MEL.m2[c(4,0)], MEL.SE2[c(4,0)])

colnames(immuno.site.df)<-c("Site","Status", "Period", "CAT", "CAT.SE", "POX", "POX.SE", "SOD", "SOD.SE", "PPO", "PPO.SE", "MEL", "MEL.SE")

Pre.immuno<-immuno.site.df[c(1:2),]
Pre.immuno$dom<-"unident"
Pre.immuno<-Pre.immuno[, c(1:3, 14, 4:13)]
Pre.immuno$int<-interaction(Pre.immuno$Site, Pre.immuno$dom)

## now this subset of site means can be overlayed with datafame from Oct 2014- Feb 2016
fig.df.immuno.comp<-rbind(Pre.immuno, fig.df.immuno)


color.scheme2<-c("gray70", 'lightskyblue', 'dodgerblue', 
                "gray40", "thistle", "coral")
break.points.immuno<-c("Lilipuna.Pre.unident", "Lilipuna.C", "Lilipuna.D", 
                "Reef 14.Pre.unident", "Reef 14.C", "Reef 14.D")
legend.names.immuno<-c("Lil.Pre unident.", "Lil C", "Lil D", 
                "R14.Pre unident.", "R14 C", "R14 D")


fig.df.immuno.comp$int<-factor(fig.df.immuno.comp$int, levels=c("Lilipuna.Pre.unident","Lilipuna.C", "Lilipuna.D", "Reef 14.Pre.unident", "Reef 14.C", "Reef 14.D"))


########## Prophenoloxidase (PPO) ############
PPO.fig<-ggplot(data=fig.df.immuno.comp, aes(x=Period, y=PPO, group=int, color=int)) + geom_errorbar(aes(ymin=PPO-PPO.SE, ymax=PPO+PPO.SE), size=.8, width=0, position=pd) +
  geom_line(position=pd, size=.8) +
  geom_point(position=pd, size=3) +
  coord_cartesian(ylim=c(0, 3)) +
  ylab(expression(paste("PO Activity", ~(Delta~Abs~"490nm"~min^-1~mg~prot^-1), sep=""))) +
  scale_color_manual(values=c(color.scheme2),
                     name= "Site:Clade",
                     breaks=break.points.immuno,
                     labels=legend.names.immuno) +
  scale_x_discrete(name ="Year and Event", 
      labels=axis.labels)+
  Fig.formatting
PPO.fig

Melanin (MEL)

MEL.lab<-aggregate(MEL~Site+Status+Period,data=data, mean)
MEL.labSE<-aggregate(MEL~Site+Status+Period,data=data, std.error)
MEL.df<-cbind(MEL.lab, MEL.labSE[c(4,0)]); colnames(MEL.df[,4:5])<-c("MEL", "MEL.labSE")
MEL.df$Period<-factor(MEL.df$Period, levels=c("2014 Field.Feb", "2014 Lab", "2014 Oct", "2015 Feb", "2015 Oct", "2016 Feb"))

plot(MEL~Period, data=MEL.df, cex.axis=0.8)

#dev.copy(pdf,"Figures/MEL.all.time.pdf", height=5, width=7, encod="MacRoman")
#dev.off

######## Melanin (MEL) #############
MEL.fig<-ggplot(data=fig.df.immuno.comp, aes(x=Period, y=MEL, group=int, color=int)) + geom_errorbar(aes(ymin=MEL-MEL.SE, ymax=MEL+MEL.SE), size=.8, width=0, position=pd) +
  geom_line(position=pd, size=.8) +
  geom_point(position=pd, size=3) +
  coord_cartesian(ylim=c(0, 0.02)) +
  ylab(expression(paste("MEL", ~(Abs~"490nm"~mg~tissue^-1), sep="")))+
  scale_color_manual(values=c(color.scheme2),
                     name= "Site:Clade",
                     breaks=break.points.immuno,
                     labels=legend.names.immuno) +
  scale_x_discrete(name ="Year and Event", 
      labels=axis.labels)+
  Fig.formatting
MEL.fig

ggsave("Figures/MEL.field.pdf", height=5, width=7, encod="MacRoman")

Peroxidase (POX)

########## Peroxidase (POX) ##############
POX.fig<-ggplot(data=fig.df.immuno.comp, aes(x=Period, y=POX, group=int, color=int)) + geom_errorbar(aes(ymin=POX-POX.SE, ymax=POX+POX.SE), size=.8, width=0, position=pd) +
  geom_line(position=pd, size=.8) +
  geom_point(position=pd, size=3) +
  coord_cartesian(ylim=c(0, 1.5)) +
  ylab(expression(paste("POX activity", ~(Delta~Abs~"470nm"~min^-1~mg~prot^-1), sep=""))) +
  scale_color_manual(values=c(color.scheme2),
                     name= "Site:Clade",
                     breaks=break.points.immuno,
                     labels=legend.names.immuno) +
  scale_x_discrete(name ="Year and Event", 
      labels=axis.labels)+
  Fig.formatting
POX.fig

Catalase (CAT)

###### Catalase (CAT) ########
CAT.fig<-ggplot(data=fig.df.immuno.comp, aes(x=Period, y=CAT, group=int, color=int)) + geom_errorbar(aes(ymin=CAT-CAT.SE, ymax=CAT+CAT.SE), size=.8, width=0, position=pd) +
  geom_line(position=pd, size=.8) +
  geom_point(position=pd, size=3) +
  coord_cartesian(ylim=c(0, 800)) +
  ylab(expression(paste("CAT activity", ~(Abs~"490nm"~mg~tissue^-1), sep="")))+
  scale_color_manual(values=c(color.scheme2),
                     name= "Site:Clade",
                     breaks=break.points.immuno,
                     labels=legend.names.immuno) +
  scale_x_discrete(name ="Year and Event", 
      labels=axis.labels)+
  Fig.formatting
CAT.fig

Superoxide dismutase (SOD)

######### Superoxide dismutase (SOD) ########
SOD.fig<-ggplot(data=fig.df.immuno.comp, aes(x=Period, y=SOD/10^3, group=int, color=int)) + geom_errorbar(aes(ymin=SOD/10^3-SOD.SE/10^3, ymax=SOD/10^3+SOD.SE/10^3), size=.8, width=0, position=pd) +
  geom_line(position=pd, size=.8) +
  geom_point(position=pd, size=3) +
  coord_cartesian(ylim=c(0, 50)) +
  ylab(expression(paste("SOD activity", ~x10^3~mg~protein^-1, sep=""))) +
  scale_color_manual(values=c(color.scheme2),
                     name= "Site:Clade",
                     breaks=break.points.immuno,
                     labels=legend.names.immuno) +
  scale_x_discrete(name ="Year and Event", 
      labels=axis.labels)+
  Fig.formatting
SOD.fig

Models Running models of response variables. First check for assumptions of ANOVA and apply transformations where necessary

df # this is the data

model.data<-df[ , !(names(df) %in% c("coral.red.adj","coral.green.adj", "coral.blue.adj", "PC1", "propC", "propD", "syms"))] # remove unwanted columns

model.data<-model.data[!c(model.data$Period=="2014 Lab"),] 
model.data<-model.data[!c(model.data$Period=="2014 Field.Feb"),] # drop data prior to Oct 2014 bleaching

model.data$Period<-factor(model.data$Period) #re-call as factor to remove dropped levels
model.data$Status<-factor(model.data$Status)


for(i in c(7:9, 11:16)){
  Y<-model.data[,i]
  full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
  R <- resid(full) #save glm residuals
  
  op<-par(mfrow = c(2,2), mar=c(5,4,1,2), pty="sq")
  plot(full, add.smooth = FALSE, which=1)
  QQ <- qqnorm(R, main = colnames(model.data)[i]) 
  QQline <- qqline(R)
  hist(R, xlab="Residuals", main = colnames(model.data)[i])
}

# need BoxCox transformation for most immuno data
Yx<-(model.data$AFDW) 

### test with transformations
boxcox(Yx~Period*Site*dom, data=model.data, lambda=seq(-0.25, 1, length=10))

### test residuals, histogram, QQ
full<-lm(Yx~Period*Site*dom, data=model.data, na.action=na.exclude)
  R <- resid(full) #save glm residuals
  op<-par(mfrow = c(2,2), mar=c(5,4,1,2), pty="sq")
  plot(full, add.smooth = FALSE, which=1)
  QQ <- qqnorm(R) 
  QQline <- qqline(R)
  hist(R, xlab="Residuals")

##### revised model.data (below are data needing transformations)
model.data$CAT<-(model.data$CAT)^0.5
model.data$POX<-(model.data$POX)^0.25
model.data$PPO<-(model.data$PPO) ^0.25
model.data$MEL<-(model.data$MEL)^0.25

## dataframe 'model.data' now ready for full model analysis

Run linear models for response variables

# Models

# symbionts ---- 
Y<-model.data$zoox
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
summary(full)
## 
## Call:
## lm(formula = Y ~ Period * Site * dom, data = model.data, na.action = na.exclude)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1991348  -350455   -68451   327150  2198987 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       634517     156105   4.065 6.14e-05 ***
## Period2015 Feb                    579684     244570   2.370 0.018408 *  
## Period2015 Oct                     71025     207209   0.343 0.732012    
## Period2016 Feb                    662557     233155   2.842 0.004794 ** 
## SiteReef 14                       193012     205163   0.941 0.347575    
## domD                             1443825     201531   7.164 6.05e-12 ***
## Period2015 Feb:SiteReef 14       -125786     313143  -0.402 0.688199    
## Period2015 Oct:SiteReef 14          2294     281470   0.008 0.993502    
## Period2016 Feb:SiteReef 14       -157264     302621  -0.520 0.603674    
## Period2015 Feb:domD             -1093240     299976  -3.644 0.000316 ***
## Period2015 Oct:domD              -422283     282317  -1.496 0.135759    
## Period2016 Feb:domD              -435738     291630  -1.494 0.136185    
## SiteReef 14:domD                 -950415     285084  -3.334 0.000964 ***
## Period2015 Feb:SiteReef 14:domD  1223395     411985   2.970 0.003223 ** 
## Period2015 Oct:SiteReef 14:domD  1101584     399311   2.759 0.006158 ** 
## Period2016 Feb:SiteReef 14:domD   655847     407091   1.611 0.108215    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 624400 on 301 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.4289, Adjusted R-squared:  0.4005 
## F-statistic: 15.07 on 15 and 301 DF,  p-value: < 2.2e-16
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                    Sum Sq  Df F value   Pr(>F)    
## Period          1.108e+13   3   9.470 5.37e-06 ***
## Site            3.068e+10   1   0.079  0.77928    
## dom             5.640e+13   1 144.647  < 2e-16 ***
## Period:Site     4.604e+12   3   3.936  0.00888 ** 
## Period:dom      3.461e+12   3   2.959  0.03260 *  
## Site:dom        8.347e+11   1   2.141  0.14447    
## Period:Site:dom 4.342e+12   3   3.712  0.01198 *  
## Residuals       1.174e+14 301                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Site*dom|Period)
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  Site     dom    lsmean       SE  df  lower.CL  upper.CL .group
##  Lilipuna C    634517.4 156105.4 301  327321.2  941713.6  a    
##  Reef 14  C    827529.1 133127.2 301  565551.3 1089506.9  ab   
##  Reef 14  D   1320939.2 151444.5 301 1022915.1 1618963.4   b   
##  Lilipuna D   2078342.6 127459.6 301 1827518.0 2329167.3    c  
## 
## Period = 2015 Feb:
##  Site     dom    lsmean       SE  df  lower.CL  upper.CL .group
##  Lilipuna C   1214201.7 188270.2 301  843709.1 1584694.3  a    
##  Reef 14  C   1281427.5 143252.2 301  999524.9 1563330.1  a    
##  Lilipuna D   1564787.2 118004.6 301 1332568.7 1797005.7  ab   
##  Reef 14  D   1904992.7 136260.0 301 1636849.9 2173135.6   b   
## 
## Period = 2015 Oct:
##  Site     dom    lsmean       SE  df  lower.CL  upper.CL .group
##  Lilipuna C    705542.2 136260.0 301  437399.4  973685.1  a    
##  Reef 14  C    900848.1 136260.0 301  632705.3 1168991.0  a    
##  Lilipuna D   1727084.3 143252.2 301 1445181.7 2008986.9   b   
##  Reef 14  D   2073558.7 143252.2 301 1791656.1 2355461.3   b   
## 
## Period = 2016 Feb:
##  Site     dom    lsmean       SE  df  lower.CL  upper.CL .group
##  Lilipuna C   1297074.8 173183.4 301  956271.2 1637878.4  a    
##  Reef 14  C   1332822.8 139624.9 301 1058058.1 1607587.4  a    
##  Reef 14  D   2046341.7 143252.2 301 1764439.1 2328244.3   b   
##  Lilipuna D   2305161.5 120170.0 301 2068681.8 2541641.3   b   
## 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="zoox", par.strip.text=list(cex=0.7))

# chla ---- 
Y<-model.data$chla
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
summary(full)
## 
## Call:
## lm(formula = Y ~ Period * Site * dom, data = model.data, na.action = na.exclude)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3721 -0.9631 -0.1189  0.9540  6.6622 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      2.68699    0.38515   6.976 1.93e-11 ***
## Period2015 Feb                   1.73762    0.60342   2.880 0.004267 ** 
## Period2015 Oct                  -0.59767    0.51124  -1.169 0.243298    
## Period2016 Feb                   1.94198    0.57525   3.376 0.000832 ***
## SiteReef 14                      0.82651    0.50619   1.633 0.103554    
## domD                             1.02348    0.49723   2.058 0.040416 *  
## Period2015 Feb:SiteReef 14      -0.67916    0.77260  -0.879 0.380074    
## Period2015 Oct:SiteReef 14      -0.17368    0.69446  -0.250 0.802682    
## Period2016 Feb:SiteReef 14      -0.19711    0.74664  -0.264 0.791968    
## Period2015 Feb:domD             -2.14866    0.74011  -2.903 0.003967 ** 
## Period2015 Oct:domD              0.27701    0.69655   0.398 0.691143    
## Period2016 Feb:domD             -1.45649    0.71952  -2.024 0.043828 *  
## SiteReef 14:domD                -0.83196    0.70337  -1.183 0.237817    
## Period2015 Feb:SiteReef 14:domD  1.95188    1.01647   1.920 0.055770 .  
## Period2015 Oct:SiteReef 14:domD  1.19736    0.98520   1.215 0.225184    
## Period2016 Feb:SiteReef 14:domD -0.01341    1.00439  -0.013 0.989355    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.541 on 301 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.2258, Adjusted R-squared:  0.1872 
## F-statistic: 5.853 on 15 and 301 DF,  p-value: 1.227e-10
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                 Sum Sq  Df F value   Pr(>F)    
## Period            85.3   3  11.986 1.96e-07 ***
## Site              23.9   1  10.049  0.00168 ** 
## dom                3.4   1   1.438  0.23137    
## Period:Site        7.0   3   0.988  0.39860    
## Period:dom        66.1   3   9.278 6.92e-06 ***
## Site:dom           0.1   1   0.029  0.86516    
## Period:Site:dom   12.7   3   1.790  0.14906    
## Residuals        714.4 301                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~dom*Period)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
##  dom Period     lsmean        SE  df lower.CL upper.CL .group
##  C   2015 Oct 2.415732 0.2377202 301 1.947928 2.883536  a    
##  C   2014 Oct 3.100247 0.2530936 301 2.602190 3.598304  ab   
##  D   2014 Oct 3.707743 0.2441870 301 3.227213 4.188273   bc  
##  D   2015 Oct 3.898917 0.2499188 301 3.407107 4.390726   bcd 
##  D   2015 Feb 3.933066 0.2223668 301 3.495476 4.370657   bcd 
##  D   2016 Feb 4.087975 0.2306646 301 3.634055 4.541894   bcd 
##  C   2015 Feb 4.498293 0.2918422 301 3.923984 5.072603    cd 
##  C   2016 Feb 4.943675 0.2744296 301 4.403631 5.483718     d 
## 
## Results are averaged over the levels of: Site 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 8 estimates 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom|Period)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  dom   lsmean        SE  df lower.CL upper.CL .group
##  C   3.100247 0.2530936 301 2.602190 3.598304  a    
##  D   3.707743 0.2441870 301 3.227213 4.188273  a    
## 
## Period = 2015 Feb:
##  dom   lsmean        SE  df lower.CL upper.CL .group
##  D   3.933066 0.2223668 301 3.495476 4.370657  a    
##  C   4.498293 0.2918422 301 3.923984 5.072603  a    
## 
## Period = 2015 Oct:
##  dom   lsmean        SE  df lower.CL upper.CL .group
##  C   2.415732 0.2377202 301 1.947928 2.883536  a    
##  D   3.898917 0.2499188 301 3.407107 4.390726   b   
## 
## Period = 2016 Feb:
##  dom   lsmean        SE  df lower.CL upper.CL .group
##  D   4.087975 0.2306646 301 3.634055 4.541894  a    
##  C   4.943675 0.2744296 301 4.403631 5.483718   b   
## 
## Results are averaged over the levels of: Site 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="chla", par.strip.text=list(cex=0.7))

# chlacell ---- 
Y<-model.data$chlacell
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
summary(full)
## 
## Call:
## lm(formula = Y ~ Period * Site * dom, data = model.data, na.action = na.exclude)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4408 -0.4285 -0.0774  0.3867  5.1924 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      4.58442    0.24173  18.965  < 2e-16 ***
## Period2015 Feb                  -0.91684    0.37872  -2.421 0.016079 *  
## Period2015 Oct                  -1.16093    0.32431  -3.580 0.000402 ***
## Period2016 Feb                  -1.16061    0.36925  -3.143 0.001840 ** 
## SiteReef 14                     -0.05916    0.31769  -0.186 0.852403    
## domD                            -2.75304    0.31207  -8.822  < 2e-16 ***
## Period2015 Feb:SiteReef 14      -0.05599    0.48490  -0.115 0.908150    
## Period2015 Oct:SiteReef 14      -0.14957    0.43840  -0.341 0.733223    
## Period2016 Feb:SiteReef 14       0.72219    0.47496   1.521 0.129437    
## Period2015 Feb:domD              1.22770    0.46451   2.643 0.008652 ** 
## Period2015 Oct:domD              1.53747    0.43971   3.497 0.000543 ***
## Period2016 Feb:domD              1.18977    0.45818   2.597 0.009878 ** 
## SiteReef 14:domD                 1.21632    0.44145   2.755 0.006225 ** 
## Period2015 Feb:SiteReef 14:domD -0.74808    0.63796  -1.173 0.241888    
## Period2015 Oct:SiteReef 14:domD -0.87229    0.62233  -1.402 0.162060    
## Period2016 Feb:SiteReef 14:domD -1.62356    0.63511  -2.556 0.011074 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9669 on 298 degrees of freedom
##   (5 observations deleted due to missingness)
## Multiple R-squared:  0.4821, Adjusted R-squared:  0.456 
## F-statistic: 18.49 on 15 and 298 DF,  p-value: < 2.2e-16
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                 Sum Sq  Df F value   Pr(>F)    
## Period           19.27   3   6.870 0.000173 ***
## Site              5.97   1   6.382 0.012044 *  
## dom             181.23   1 193.843  < 2e-16 ***
## Period:Site       4.06   3   1.449 0.228640    
## Period:dom       14.11   3   5.029 0.002051 ** 
## Site:dom          3.23   1   3.459 0.063903 .  
## Period:Site:dom   6.15   3   2.192 0.089017 .  
## Residuals       278.61 298                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~dom*Period)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
##  dom Period     lsmean        SE  df lower.CL upper.CL .group
##  D   2016 Feb 1.988445 0.1447701 298 1.703544 2.273347  a    
##  D   2015 Oct 2.275578 0.1590180 298 1.962637 2.588518  a    
##  D   2015 Feb 2.318790 0.1395622 298 2.044137 2.593442  a    
##  D   2014 Oct 2.409966 0.1532570 298 2.108363 2.711569  a    
##  C   2015 Oct 3.319127 0.1510518 298 3.021863 3.616390   b   
##  C   2015 Feb 3.610003 0.1831665 298 3.249540 3.970467   b   
##  C   2016 Feb 3.755331 0.1765338 298 3.407920 4.102742   b   
##  C   2014 Oct 4.554842 0.1588470 298 4.242238 4.867446    c  
## 
## Results are averaged over the levels of: Site 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 8 estimates 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom|Period)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  dom   lsmean        SE  df lower.CL upper.CL .group
##  D   2.409966 0.1532570 298 2.108363 2.711569  a    
##  C   4.554842 0.1588470 298 4.242238 4.867446   b   
## 
## Period = 2015 Feb:
##  dom   lsmean        SE  df lower.CL upper.CL .group
##  D   2.318790 0.1395622 298 2.044137 2.593442  a    
##  C   3.610003 0.1831665 298 3.249540 3.970467   b   
## 
## Period = 2015 Oct:
##  dom   lsmean        SE  df lower.CL upper.CL .group
##  D   2.275578 0.1590180 298 1.962637 2.588518  a    
##  C   3.319127 0.1510518 298 3.021863 3.616390   b   
## 
## Period = 2016 Feb:
##  dom   lsmean        SE  df lower.CL upper.CL .group
##  D   1.988445 0.1447701 298 1.703544 2.273347  a    
##  C   3.755331 0.1765338 298 3.407920 4.102742   b   
## 
## Results are averaged over the levels of: Site 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="chlacell", par.strip.text=list(cex=0.7))

# AFDW ---- 
Y<-model.data$AFDW
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
summary(full)
## 
## Call:
## lm(formula = Y ~ Period * Site * dom, data = model.data, na.action = na.exclude)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -21.358  -4.662  -0.560   4.090  47.519 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      22.5410     2.0286  11.112  < 2e-16 ***
## Period2015 Feb                   -9.8766     3.1782  -3.108  0.00207 ** 
## Period2015 Oct                    6.9224     2.6927   2.571  0.01063 *  
## Period2016 Feb                    6.3194     3.0299   2.086  0.03785 *  
## SiteReef 14                      -3.5154     2.6661  -1.319  0.18832    
## domD                              3.6565     2.6189   1.396  0.16369    
## Period2015 Feb:SiteReef 14        7.6903     4.0693   1.890  0.05974 .  
## Period2015 Oct:SiteReef 14        0.9363     3.6577   0.256  0.79814    
## Period2016 Feb:SiteReef 14       -5.2741     3.9326  -1.341  0.18089    
## Period2015 Feb:domD              -2.8045     3.8982  -0.719  0.47242    
## Period2015 Oct:domD              11.6543     3.6687   3.177  0.00164 ** 
## Period2016 Feb:domD              -2.4168     3.7898  -0.638  0.52414    
## SiteReef 14:domD                 -4.9778     3.7047  -1.344  0.18007    
## Period2015 Feb:SiteReef 14:domD   4.5821     5.3538   0.856  0.39275    
## Period2015 Oct:SiteReef 14:domD  -4.5805     5.1891  -0.883  0.37809    
## Period2016 Feb:SiteReef 14:domD   4.8994     5.2902   0.926  0.35512    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.114 on 301 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.5017, Adjusted R-squared:  0.4769 
## F-statistic:  20.2 on 15 and 301 DF,  p-value: < 2.2e-16
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                 Sum Sq  Df F value   Pr(>F)    
## Period           14064   3  71.198  < 2e-16 ***
## Site              1616   1  24.536 1.22e-06 ***
## dom                937   1  14.228 0.000195 ***
## Period:Site       1999   3  10.118 2.27e-06 ***
## Period:dom        1303   3   6.599 0.000248 ***
## Site:dom           288   1   4.381 0.037182 *  
## Period:Site:dom    287   3   1.452 0.227787    
## Residuals        19819 301                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Site*dom)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
##  Site     dom   lsmean        SE  df lower.CL upper.CL .group
##  Reef 14  C   20.70499 0.8974235 301 18.93897 22.47101  a    
##  Reef 14  D   22.21717 0.9333888 301 20.38037 24.05396  a    
##  Lilipuna C   23.38230 1.0694994 301 21.27765 25.48694  a    
##  Lilipuna D   28.64702 0.8291230 301 27.01541 30.27863   b   
## 
## Results are averaged over the levels of: Period 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site|Period)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  Site       lsmean       SE  df lower.CL upper.CL .group
##  Reef 14  18.36490 1.310155 301 15.78667 20.94312  a    
##  Lilipuna 24.36925 1.309455 301 21.79240 26.94609   b   
## 
## Period = 2015 Feb:
##  Site       lsmean       SE  df lower.CL upper.CL .group
##  Lilipuna 13.09036 1.443721 301 10.24930 15.93143  a    
##  Reef 14  17.06741 1.284606 301 14.53946 19.59535   b   
## 
## Period = 2015 Oct:
##  Site       lsmean       SE  df lower.CL upper.CL .group
##  Reef 14  29.76049 1.284606 301 27.23254 32.28843  a    
##  Lilipuna 37.11879 1.284606 301 34.59085 39.64674   b   
## 
## Period = 2016 Feb:
##  Site       lsmean       SE  df lower.CL upper.CL .group
##  Reef 14  20.65152 1.299771 301 18.09373 23.20931  a    
##  Lilipuna 29.48023 1.369628 301 26.78497 32.17549   b   
## 
## Results are averaged over the levels of: dom 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom|Period)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  dom   lsmean       SE  df lower.CL upper.CL .group
##  C   20.78329 1.333051 301 18.16001 23.40657  a    
##  D   21.95086 1.286140 301 19.41989 24.48182  a    
## 
## Period = 2015 Feb:
##  dom   lsmean       SE  df lower.CL upper.CL .group
##  C   14.75184 1.537141 301 11.72694 17.77674  a    
##  D   15.40593 1.171212 301 13.10113 17.71073  a    
## 
## Period = 2015 Oct:
##  dom   lsmean       SE  df lower.CL upper.CL .group
##  C   28.17382 1.252079 301 25.70988 30.63775  a    
##  D   38.70546 1.316329 301 36.11509 41.29584   b   
## 
## Period = 2016 Feb:
##  dom   lsmean       SE  df lower.CL upper.CL .group
##  C   24.46563 1.445428 301 21.62121 27.31006  a    
##  D   25.66612 1.214917 301 23.27531 28.05692  a    
## 
## Results are averaged over the levels of: Site 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="AFDW", par.strip.text=list(cex=0.7))

# prot ---- 
Y<-model.data$prot
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
summary(full)
## 
## Call:
## lm(formula = Y ~ Period * Site * dom, data = model.data, na.action = na.exclude)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.34266 -0.11312 -0.02087  0.09214  0.54100 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      0.49861    0.04088  12.198   <2e-16 ***
## Period2015 Feb                  -0.09029    0.06404  -1.410   0.1596    
## Period2015 Oct                  -0.08976    0.05426  -1.654   0.0991 .  
## Period2016 Feb                  -0.01735    0.06105  -0.284   0.7764    
## SiteReef 14                     -0.08663    0.05372  -1.613   0.1079    
## domD                             0.08155    0.05277   1.545   0.1233    
## Period2015 Feb:SiteReef 14       0.15566    0.08200   1.898   0.0586 .  
## Period2015 Oct:SiteReef 14       0.08961    0.07370   1.216   0.2250    
## Period2016 Feb:SiteReef 14       0.03516    0.07924   0.444   0.6576    
## Period2015 Feb:domD             -0.07079    0.07855  -0.901   0.3682    
## Period2015 Oct:domD              0.01911    0.07393   0.259   0.7962    
## Period2016 Feb:domD             -0.07189    0.07636  -0.941   0.3473    
## SiteReef 14:domD                -0.08884    0.07531  -1.180   0.2390    
## Period2015 Feb:SiteReef 14:domD  0.09762    0.10833   0.901   0.3683    
## Period2015 Oct:SiteReef 14:domD  0.03793    0.10503   0.361   0.7183    
## Period2016 Feb:SiteReef 14:domD  0.13003    0.10706   1.215   0.2255    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1635 on 300 degrees of freedom
##   (3 observations deleted due to missingness)
## Multiple R-squared:  0.08939,    Adjusted R-squared:  0.04385 
## F-statistic: 1.963 on 15 and 300 DF,  p-value: 0.01762
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                 Sum Sq  Df F value  Pr(>F)   
## Period           0.058   3   0.722 0.53971   
## Site             0.053   1   1.980 0.16042   
## dom              0.126   1   4.706 0.03085 * 
## Period:Site      0.412   3   5.135 0.00178 **
## Period:dom       0.035   3   0.440 0.72477   
## Site:dom         0.011   1   0.394 0.53082   
## Period:Site:dom  0.048   3   0.596 0.61820   
## Residuals        8.020 300                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~dom)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
##  dom    lsmean         SE  df  lower.CL  upper.CL .group
##  C   0.4410006 0.01406619 300 0.4133197 0.4686815  a    
##  D   0.4804365 0.01263926 300 0.4555637 0.5053093   b   
## 
## Results are averaged over the levels of: Period, Site 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site|Period)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  Site        lsmean         SE  df  lower.CL  upper.CL .group
##  Reef 14  0.4083375 0.02686119 300 0.3554773 0.4611977  a    
##  Lilipuna 0.5393909 0.02638573 300 0.4874664 0.5913155   b   
## 
## Period = 2015 Feb:
##  Site        lsmean         SE  df  lower.CL  upper.CL .group
##  Lilipuna 0.4137070 0.02909120 300 0.3564583 0.4709556  a    
##  Reef 14  0.4871237 0.02588501 300 0.4361845 0.5380629  a    
## 
## Period = 2015 Oct:
##  Site        lsmean         SE  df  lower.CL  upper.CL .group
##  Reef 14  0.4367003 0.02588501 300 0.3857612 0.4876395  a    
##  Lilipuna 0.4591846 0.02588501 300 0.4082454 0.5101238  a    
## 
## Period = 2016 Feb:
##  Site        lsmean         SE  df  lower.CL  upper.CL .group
##  Reef 14  0.4552101 0.02619059 300 0.4036696 0.5067507  a    
##  Lilipuna 0.4860943 0.02759823 300 0.4317837 0.5404049  a    
## 
## Results are averaged over the levels of: dom 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="prot", par.strip.text=list(cex=0.7))

######################## 
### immunology
########################

# CAT ---- 
Y<-model.data$CAT
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
summary(full)
## 
## Call:
## lm(formula = Y ~ Period * Site * dom, data = model.data, na.action = na.exclude)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.1768  -1.7001   0.1203   1.7657   9.5328 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      11.5496     0.7366  15.679  < 2e-16 ***
## Period2015 Feb                   -3.1503     1.1541  -2.730  0.00672 ** 
## Period2015 Oct                    7.4748     0.9883   7.563 5.17e-13 ***
## Period2016 Feb                   -1.2333     1.1002  -1.121  0.26322    
## SiteReef 14                      -0.1867     0.9778  -0.191  0.84870    
## domD                              2.7246     0.9592   2.840  0.00482 ** 
## Period2015 Feb:SiteReef 14        0.1926     1.4840   0.130  0.89682    
## Period2015 Oct:SiteReef 14       -3.6517     1.3590  -2.687  0.00763 ** 
## Period2016 Feb:SiteReef 14       -1.8530     1.4345  -1.292  0.19749    
## Period2015 Feb:domD              -1.1980     1.4210  -0.843  0.39989    
## Period2015 Oct:domD              -1.9804     1.3656  -1.450  0.14808    
## Period2016 Feb:domD              -3.0389     1.3818  -2.199  0.02865 *  
## SiteReef 14:domD                 -1.5819     1.3580  -1.165  0.24504    
## Period2015 Feb:SiteReef 14:domD   2.3155     1.9529   1.186  0.23672    
## Period2015 Oct:SiteReef 14:domD   0.4003     1.9464   0.206  0.83720    
## Period2016 Feb:SiteReef 14:domD   2.4720     1.9299   1.281  0.20125    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.946 on 291 degrees of freedom
##   (12 observations deleted due to missingness)
## Multiple R-squared:  0.5888, Adjusted R-squared:  0.5676 
## F-statistic: 27.78 on 15 and 291 DF,  p-value: < 2.2e-16
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                 Sum Sq  Df F value   Pr(>F)    
## Period          3048.5   3 117.049  < 2e-16 ***
## Site             185.2   1  21.334 5.80e-06 ***
## dom               81.0   1   9.328  0.00247 ** 
## Period:Site      228.3   3   8.766 1.39e-05 ***
## Period:dom        57.9   3   2.222  0.08568 .  
## Site:dom           1.6   1   0.183  0.66900    
## Period:Site:dom   22.5   3   0.862  0.46111    
## Residuals       2526.3 291                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~dom)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
##  dom   lsmean        SE  df lower.CL upper.CL .group
##  C   11.56504 0.2560231 291 11.06115 12.06893  a    
##  D   12.59285 0.2319532 291 12.13633 13.04937   b   
## 
## Results are averaged over the levels of: Period, Site 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Period)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
##  Period      lsmean        SE  df  lower.CL upper.CL .group
##  2015 Feb  9.348931 0.3508581 291  8.658390 10.03947  a    
##  2016 Feb  9.361876 0.3428152 291  8.687164 10.03659  a    
##  2014 Oct 12.423079 0.3394984 291 11.754895 13.09126   b   
##  2015 Oct 17.181893 0.3485938 291 16.495809 17.86798    c  
## 
## Results are averaged over the levels of: Site, dom 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site|Period)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  Site        lsmean        SE  df  lower.CL  upper.CL .group
##  Reef 14  11.934266 0.4806477 291 10.988279 12.880252  a    
##  Lilipuna 12.911891 0.4795982 291 11.967970 13.855812  a    
## 
## Period = 2015 Feb:
##  Site        lsmean        SE  df  lower.CL  upper.CL .group
##  Lilipuna  9.162569 0.5242351 291  8.130796 10.194342  a    
##  Reef 14   9.535293 0.4664582 291  8.617233 10.453352  a    
## 
## Period = 2015 Oct:
##  Site        lsmean        SE  df  lower.CL  upper.CL .group
##  Reef 14  14.967313 0.4998807 291 13.983473 15.951153  a    
##  Lilipuna 19.396474 0.4859936 291 18.439966 20.352982   b   
## 
## Period = 2016 Feb:
##  Site        lsmean        SE  df  lower.CL  upper.CL .group
##  Reef 14   8.564590 0.4719648 291  7.635693  9.493487  a    
##  Lilipuna 10.159162 0.4973311 291  9.180340 11.137984   b   
## 
## Results are averaged over the levels of: dom 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="CAT", par.strip.text=list(cex=0.7))

# POX ---- 
Y<-model.data$POX
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
summary(full)
## 
## Call:
## lm(formula = Y ~ Period * Site * dom, data = model.data, na.action = na.exclude)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.65920 -0.11092 -0.00442  0.10684  0.58141 
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      0.8592630  0.0472279  18.194  < 2e-16 ***
## Period2015 Feb                   0.0008085  0.0726087   0.011  0.99112    
## Period2015 Oct                  -0.1868710  0.0647961  -2.884  0.00423 ** 
## Period2016 Feb                   0.0951206  0.0693115   1.372  0.17103    
## SiteReef 14                      0.0137554  0.0612474   0.225  0.82246    
## domD                            -0.0352386  0.0602039  -0.585  0.55880    
## Period2015 Feb:SiteReef 14      -0.0243207  0.0924862  -0.263  0.79277    
## Period2015 Oct:SiteReef 14      -0.0750079  0.0864882  -0.867  0.38653    
## Period2016 Feb:SiteReef 14      -0.0557302  0.0894300  -0.623  0.53367    
## Period2015 Feb:domD              0.0038584  0.0886621   0.044  0.96532    
## Period2015 Oct:domD              0.1738882  0.0876564   1.984  0.04825 *  
## Period2016 Feb:domD              0.0200563  0.0862397   0.233  0.81627    
## SiteReef 14:domD                -0.0529785  0.0858822  -0.617  0.53781    
## Period2015 Feb:SiteReef 14:domD  0.0669825  0.1226620   0.546  0.58544    
## Period2015 Oct:SiteReef 14:domD  0.0213673  0.1242025   0.172  0.86353    
## Period2016 Feb:SiteReef 14:domD -0.0029081  0.1209226  -0.024  0.98083    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1829 on 284 degrees of freedom
##   (19 observations deleted due to missingness)
## Multiple R-squared:  0.1939, Adjusted R-squared:  0.1513 
## F-statistic: 4.554 on 15 and 284 DF,  p-value: 9.044e-08
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                 Sum Sq  Df F value   Pr(>F)    
## Period           1.657   3  16.504 6.51e-10 ***
## Site             0.121   1   3.619   0.0581 .  
## dom              0.001   1   0.042   0.8382    
## Period:Site      0.089   3   0.884   0.4500    
## Period:dom       0.363   3   3.612   0.0138 *  
## Site:dom         0.018   1   0.547   0.4602    
## Period:Site:dom  0.014   3   0.138   0.9372    
## Residuals        9.502 284                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Period)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
##  Period      lsmean         SE  df  lower.CL  upper.CL .group
##  2015 Oct 0.7031878 0.02243114 284 0.6590354 0.7473401  a    
##  2014 Oct 0.8352768 0.02147056 284 0.7930152 0.8775384   b   
##  2015 Feb 0.8425998 0.02189491 284 0.7995029 0.8856967   b   
##  2016 Feb 0.9118334 0.02128162 284 0.8699437 0.9537232   b   
## 
## Results are averaged over the levels of: Site, dom 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom|Period)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  dom    lsmean         SE  df  lower.CL  upper.CL .group
##  D   0.8044129 0.03010197 284 0.7451616 0.8636641  a    
##  C   0.8661407 0.03062371 284 0.8058625 0.9264190  a    
## 
## Period = 2015 Feb:
##  dom    lsmean         SE  df  lower.CL  upper.CL .group
##  D   0.8304106 0.02677568 284 0.7777067 0.8831146  a    
##  C   0.8547889 0.03464985 284 0.7865858 0.9229920  a    
## 
## Period = 2015 Oct:
##  dom    lsmean         SE  df  lower.CL  upper.CL .group
##  C   0.6417658 0.03053261 284 0.5816668 0.7018647  a    
##  D   0.7646097 0.03286920 284 0.6999116 0.8293079   b   
## 
## Period = 2016 Feb:
##  dom    lsmean         SE  df  lower.CL  upper.CL .group
##  D   0.8902706 0.02738635 284 0.8363646 0.9441766  a    
##  C   0.9333963 0.03258248 284 0.8692625 0.9975301  a    
## 
## Results are averaged over the levels of: Site 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site|Period)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  Site        lsmean         SE  df  lower.CL  upper.CL .group
##  Reef 14  0.8289099 0.03062371 284 0.7686316 0.8891881  a    
##  Lilipuna 0.8416437 0.03010197 284 0.7823925 0.9008950  a    
## 
## Period = 2015 Feb:
##  Site        lsmean         SE  df  lower.CL  upper.CL .group
##  Reef 14  0.8408181 0.02929911 284 0.7831471 0.8984890  a    
##  Lilipuna 0.8443814 0.03254399 284 0.7803234 0.9084395  a    
## 
## Period = 2015 Oct:
##  Site        lsmean         SE  df  lower.CL  upper.CL .group
##  Reef 14  0.6646587 0.03158864 284 0.6024811 0.7268362  a    
##  Lilipuna 0.7417168 0.03185564 284 0.6790137 0.8044200  a    
## 
## Period = 2016 Feb:
##  Site        lsmean         SE  df  lower.CL  upper.CL .group
##  Reef 14  0.8768744 0.02929911 284 0.8192034 0.9345453  a    
##  Lilipuna 0.9467925 0.03087382 284 0.8860219 1.0075630  a    
## 
## Results are averaged over the levels of: dom 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="POX", par.strip.text=list(cex=0.7))

# SOD ---- 
Y<-model.data$SOD
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
summary(full)
## 
## Call:
## lm(formula = Y ~ Period * Site * dom, data = model.data, na.action = na.exclude)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -20899.7  -4928.7   -938.1   4206.2  26108.1 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      17047.3     1838.1   9.275  < 2e-16 ***
## Period2015 Feb                    4690.1     2879.7   1.629   0.1044    
## Period2015 Oct                   13071.0     2439.8   5.357 1.69e-07 ***
## Period2016 Feb                   16540.8     2745.3   6.025 4.96e-09 ***
## SiteReef 14                      -3050.5     2415.7  -1.263   0.2076    
## domD                             -1452.9     2372.9  -0.612   0.5408    
## Period2015 Feb:SiteReef 14         700.1     3687.1   0.190   0.8495    
## Period2015 Oct:SiteReef 14        2107.8     3314.2   0.636   0.5253    
## Period2016 Feb:SiteReef 14        -191.6     3563.2  -0.054   0.9572    
## Period2015 Feb:domD               6526.6     3532.1   1.848   0.0656 .  
## Period2015 Oct:domD                924.7     3347.8   0.276   0.7826    
## Period2016 Feb:domD                738.9     3433.8   0.215   0.8298    
## SiteReef 14:domD                  2621.0     3356.7   0.781   0.4355    
## Period2015 Feb:SiteReef 14:domD  -4670.4     4850.9  -0.963   0.3364    
## Period2015 Oct:SiteReef 14:domD  -3039.1     4735.2  -0.642   0.5215    
## Period2016 Feb:SiteReef 14:domD   1278.3     4793.3   0.267   0.7899    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7352 on 299 degrees of freedom
##   (4 observations deleted due to missingness)
## Multiple R-squared:  0.4722, Adjusted R-squared:  0.4457 
## F-statistic: 17.83 on 15 and 299 DF,  p-value: < 2.2e-16
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                    Sum Sq  Df F value Pr(>F)    
## Period          1.349e+10   3  83.207 <2e-16 ***
## Site            2.631e+08   1   4.867 0.0281 *  
## dom             8.110e+07   1   1.500 0.2216    
## Period:Site     9.381e+07   3   0.578 0.6296    
## Period:dom      2.319e+08   3   1.430 0.2340    
## Site:dom        2.041e+07   1   0.378 0.5394    
## Period:Site:dom 1.021e+08   3   0.630 0.5964    
## Residuals       1.616e+10 299                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Period)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
##  Period     lsmean       SE  df lower.CL upper.CL .group
##  2014 Oct 15450.86 839.1840 299 13799.40 17102.31  a    
##  2015 Feb 22586.64 875.4958 299 20863.73 24309.56   b   
##  2015 Oct 29278.28 834.9551 299 27635.15 30921.42    c  
##  2016 Feb 32584.88 855.4263 299 30901.46 34268.30     d 
## 
## Results are averaged over the levels of: Site, dom 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
##  Site       lsmean       SE  df lower.CL upper.CL .group
##  Reef 14  24030.26 588.7094 299 22871.72 25188.80  a    
##  Lilipuna 25920.07 615.0838 299 24709.63 27130.52   b   
## 
## Results are averaged over the levels of: Period, dom 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site|Period)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  Site       lsmean       SE  df lower.CL upper.CL .group
##  Reef 14  14580.86 1187.103 299 12244.73 16917.00  a    
##  Lilipuna 16320.86 1186.468 299 13985.97 18655.74  a    
## 
## Period = 2015 Feb:
##  Site       lsmean       SE  df lower.CL upper.CL .group
##  Reef 14  20899.09 1163.953 299 18608.52 23189.67  a    
##  Lilipuna 24274.19 1308.123 299 21699.90 26848.49  a    
## 
## Period = 2015 Oct:
##  Site       lsmean       SE  df lower.CL upper.CL .group
##  Reef 14  28702.40 1180.805 299 26378.66 31026.14  a    
##  Lilipuna 29854.17 1180.805 299 27530.43 32177.91  a    
## 
## Period = 2016 Feb:
##  Site       lsmean       SE  df lower.CL upper.CL .group
##  Reef 14  31938.68 1177.693 299 29621.06 34256.30  a    
##  Lilipuna 33231.08 1240.990 299 30788.90 35673.26  a    
## 
## Results are averaged over the levels of: dom 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="SOD", par.strip.text=list(cex=0.7))

# PPO ---- 
Y<-model.data$PPO
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
summary(full)
## 
## Call:
## lm(formula = Y ~ Period * Site * dom, data = model.data, na.action = na.exclude)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.26790 -0.07227 -0.01267  0.06667  0.48933 
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      0.7029273  0.0294739  23.849  < 2e-16 ***
## Period2015 Feb                   0.3324377  0.0453135   7.336 2.12e-12 ***
## Period2015 Oct                  -0.0715260  0.0385904  -1.853   0.0648 .  
## Period2016 Feb                   0.0229863  0.0432558   0.531   0.5955    
## SiteReef 14                      0.0346904  0.0382232   0.908   0.3648    
## domD                             0.0053348  0.0375720   0.142   0.8872    
## Period2015 Feb:SiteReef 14      -0.0008834  0.0577186  -0.015   0.9878    
## Period2015 Oct:SiteReef 14      -0.0099309  0.0519810  -0.191   0.8486    
## Period2016 Feb:SiteReef 14      -0.0151842  0.0558113  -0.272   0.7858    
## Period2015 Feb:domD              0.0413138  0.0553321   0.747   0.4559    
## Period2015 Oct:domD              0.0139221  0.0524984   0.265   0.7910    
## Period2016 Feb:domD              0.0370833  0.0541778   0.684   0.4942    
## SiteReef 14:domD                 0.0018843  0.0526352   0.036   0.9715    
## Period2015 Feb:SiteReef 14:domD -0.0155636  0.0756755  -0.206   0.8372    
## Period2015 Oct:SiteReef 14:domD -0.0065871  0.0738872  -0.089   0.9290    
## Period2016 Feb:SiteReef 14:domD  0.0065517  0.0750427   0.087   0.9305    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1142 on 296 degrees of freedom
##   (7 observations deleted due to missingness)
## Multiple R-squared:  0.6877, Adjusted R-squared:  0.6719 
## F-statistic: 43.45 on 15 and 296 DF,  p-value: < 2.2e-16
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                 Sum Sq  Df F value Pr(>F)    
## Period           8.112   3 207.503 <2e-16 ***
## Site             0.055   1   4.227 0.0407 *  
## dom              0.054   1   4.135 0.0429 *  
## Period:Site      0.002   3   0.051 0.9848    
## Period:dom       0.020   3   0.510 0.6757    
## Site:dom         0.000   1   0.005 0.9419    
## Period:Site:dom  0.001   3   0.031 0.9927    
## Residuals        3.857 296                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~dom)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
##  dom    lsmean          SE  df  lower.CL  upper.CL .group
##  C   0.7879973 0.009863396 296 0.7685860 0.8074085  a    
##  D   0.8154042 0.008883127 296 0.7979221 0.8328863   b   
## 
## Results are averaged over the levels of: Period, Site 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Period)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
##  Period      lsmean         SE  df  lower.CL  upper.CL .group
##  2015 Oct 0.6522339 0.01296356 296 0.6267215 0.6777463  a    
##  2014 Oct 0.7234110 0.01315880 296 0.6975144 0.7493077   b   
##  2016 Feb 0.7589848 0.01337193 296 0.7326687 0.7853009   b   
##  2015 Feb 1.0721731 0.01359300 296 1.0454219 1.0989242    c  
## 
## Results are averaged over the levels of: Site, dom 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="PPO", par.strip.text=list(cex=0.7))

# MEL ---- 
Y<-model.data$MEL
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
summary(full)
## 
## Call:
## lm(formula = Y ~ Period * Site * dom, data = model.data, na.action = na.exclude)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.085284 -0.014451 -0.000336  0.011635  0.188528 
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      0.3454562  0.0056104  61.574   <2e-16 ***
## Period2015 Feb                  -0.2169668  0.0087898 -24.684   <2e-16 ***
## Period2015 Oct                  -0.1051567  0.0075272 -13.970   <2e-16 ***
## Period2016 Feb                  -0.1352463  0.0083796 -16.140   <2e-16 ***
## SiteReef 14                     -0.0033742  0.0073735  -0.458   0.6476    
## domD                            -0.0001185  0.0073057  -0.016   0.9871    
## Period2015 Feb:SiteReef 14       0.0279336  0.0112543   2.482   0.0136 *  
## Period2015 Oct:SiteReef 14      -0.0152304  0.0101751  -1.497   0.1355    
## Period2016 Feb:SiteReef 14       0.0073611  0.0108762   0.677   0.4991    
## Period2015 Feb:domD             -0.0052271  0.0108233  -0.483   0.6295    
## Period2015 Oct:domD             -0.0022054  0.0103216  -0.214   0.8310    
## Period2016 Feb:domD              0.0033341  0.0105246   0.317   0.7516    
## SiteReef 14:domD                 0.0035193  0.0102903   0.342   0.7326    
## Period2015 Feb:SiteReef 14:domD -0.0071391  0.0148375  -0.481   0.6308    
## Period2015 Oct:SiteReef 14:domD -0.0012880  0.0145263  -0.089   0.9294    
## Period2016 Feb:SiteReef 14:domD  0.0033165  0.0146620   0.226   0.8212    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02244 on 297 degrees of freedom
##   (6 observations deleted due to missingness)
## Multiple R-squared:  0.9211, Adjusted R-squared:  0.9171 
## F-statistic: 231.2 on 15 and 297 DF,  p-value: < 2.2e-16
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
## 
## Response: Y
##                 Sum Sq  Df  F value   Pr(>F)    
## Period          1.7128   3 1133.636  < 2e-16 ***
## Site            0.0006   1    1.112    0.292    
## dom             0.0000   1    0.000    0.987    
## Period:Site     0.0155   3   10.241 1.95e-06 ***
## Period:dom      0.0019   3    1.276    0.283    
## Site:dom        0.0001   1    0.198    0.657    
## Period:Site:dom 0.0003   3    0.170    0.917    
## Residuals       0.1496 297                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Period)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
##  Period      lsmean          SE  df  lower.CL  upper.CL .group
##  2015 Feb 0.1371913 0.002672314 297 0.1319323 0.1424504  a    
##  2016 Feb 0.2155202 0.002611055 297 0.2103817 0.2206587   b   
##  2015 Oct 0.2303931 0.002563231 297 0.2253487 0.2354375    c  
##  2014 Oct 0.3445897 0.002572585 297 0.3395269 0.3496525     d 
## 
## Results are averaged over the levels of: Site, dom 
## Confidence level used: 0.95 
## P value adjustment: tukey method for comparing a family of 4 estimates 
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site|Period)
## NOTE: Results may be misleading due to involvement in interactions
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
##  Site        lsmean          SE  df  lower.CL  upper.CL .group
##  Reef 14  0.3437824 0.003623445 297 0.3366515 0.3509133  a    
##  Lilipuna 0.3453970 0.003652864 297 0.3382082 0.3525857  a    
## 
## Period = 2015 Feb:
##  Site        lsmean          SE  df  lower.CL  upper.CL .group
##  Lilipuna 0.1258166 0.003992841 297 0.1179588 0.1336745  a    
##  Reef 14  0.1485661 0.003552783 297 0.1415743 0.1555579   b   
## 
## Period = 2015 Oct:
##  Site        lsmean          SE  df  lower.CL  upper.CL .group
##  Reef 14  0.2216486 0.003604222 297 0.2145556 0.2287417  a    
##  Lilipuna 0.2391376 0.003645572 297 0.2319632 0.2463120   b   
## 
## Period = 2016 Feb:
##  Site        lsmean          SE  df  lower.CL  upper.CL .group
##  Lilipuna 0.2118178 0.003787927 297 0.2043632 0.2192723  a    
##  Reef 14  0.2192226 0.003594725 297 0.2121482 0.2262969  a    
## 
## Results are averaged over the levels of: dom 
## Confidence level used: 0.95 
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="prot", par.strip.text=list(cex=0.7))

### symbiodinium data


###############################
## making figures, tables, analyses

###############################
# make a table by dominant symb

df

str(df)
print(levels(df$Period))

# 4 symbiont categories
symbcomp=table(df$syms, df$Period:df$Site) 
prop.table(symbcomp, margin = 2)

# 2 symbiont categories
domsymb=table(df$dom, df$Period:df$Site)
prop.table(domsymb, margin = 2)


# to specify order in figure, can use... ...(prop.table(symbcomp[,c(1,3,2,4)], 3,4,1,2,5,6
barplot(prop.table(symbcomp, margin = 2), col = c("slategray4", "slategray2", "dodgerblue", "darkturquoise"), main= "2014 - 2016 qPCR", xlab = "Site and Status", cex.names=0.65, ylab = "Proportion of Symbiont Types")
legend("topright", cex=0.8, bty="n", legend=c("C","CD", "D", "DC"), fill=c("slategray4", "slategray2", "dodgerblue", "darkturquoise")) #inset = c(-.25, 0), xpd = NA, x.intersp=0.1, y.intersp=0.7)

dev.copy(pdf, "Figures/symb_4 levels.pdf", encod="MacRoman", height=6, width=14)
dev.off()


#####################
## figures for 2014-2016 dominant symbiont composition figure (2 categories)

# to sepcify order in figure, can use... ...(prop.table(symbcomp[,c(1,3,2,4)],
barplot(prop.table(domsymb, margin = 2), col = c("slategray4", "darkturquoise"), main= "2014 - 2015 qPCR", xlab = "Site and Status", cex.names=0.65, ylab = "Proportion of Symbiont Types")
legend("topright", legend=c("C", "D"), bty="n", fill=c("slategray4", "darkturquoise"), inset = c(-.23, 0), xpd = NA, x.intersp=0.1, y.intersp=0.7)

dev.copy(pdf, "Figures/symb_2 levels.pdf", encod="MacRoman", height=6, width=14)
dev.off()

#Chi Squared test for independence
chisq.test(symbcomp) # p<0.01, accept Ha that symcomb IS related to events
chisq.test(domsymb) # p=0.1, accept Ho that domsymb is NOT related to events

prop.table(symbcomp, margin = 2)
par(mar=c(3, 4, 2, 6))


# to specify order in figure, can use... ...(prop.table(symbcomp[,c(1,3,2,4)],
barplot(prop.table(symbcomp, margin = 2), col = c("slategray4", "slategray2", "dodgerblue", "darkturquoise"), main= "2015 - 2016 qPCR", xlab = "Event and Site", cex.names=0.65, ylab = "Proportion of Symbiont Types")
legend("topright", legend=c("C","CD", "D", "DC"), fill=c("slategray4", "slategray2", "dodgerblue", "darkturquoise"))
#inset = c(-.15, 0), xpd = NA)

##########
prop.table(domsymb, margin = 2)
par(mar=c(3, 4, 2, 6))

barplot(prop.table(domsymb, margin = 2), col = c("slategray4", "darkturquoise"), main= "2015 - 2016 qPCR", xlab = "Site and Status", cex.names=0.65, ylab = "Proportion of Symbiont Types")
legend("topright", legend=c("C", "D"), fill=c("slategray4", "darkturquoise"), inset = c(-.15, 0), xpd = NA)